



PINES* NONSIHGULAR GRAVITATIONAL POTENTIAL 
DERIVATION, DESCRIPTION AND IMPLEMENTATION 


9 FEBRUARY 1976 


HDC WOOl 3 


JOSE L. SPENCER 


MCDONNELL DOUGLAS TECHNICAL SERVICES COMPANY, INC. 
HOUSTON ASTRONAUTICS DIVISION 
HOUSTON, TEXAS 77058 (713 488-5660) 


PREPARED FOR NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 
LYNDON B. JOHNSON SPACE CENTER 
HOUSTON, TEXAS 77058 


CONTRACT NAS 9-13970 



This report was prepared by the McDonnell Douglas Technical 
Services Connany, Inc., under contract NAS 9-13970, Snace 
Shuttle Enaineerina and Operations Support, for the Johnson 
Space Center of the National Aeronautics and Snace Adminis- 
tration, Reginald M. Machell , Technical Manager. This work 
was administered under the technical direction of the Mathe- 
matical Physics Branch, Mission Planning and Analysis Division, 
Data Systems and Analysis Directorate v/ith James C. Kirkpatrick 
as Technical Monitor. 



MDC W001 3 
9 February 1976 


CONTENTS 

Section Page 

1.0 SUMMARY 1 

2-0 INTRODUCTION 1 

3.0 DISCUSSION 2 

3.1 Statement of the Problem 2 

3*2 The Proposed Method of Solution 5 

3*3 The "Derived” Legendre Functions 7 

3*4 The Complex Representation 3 

3-5 The Recurrence Process 10 

. 3*5.1 Complex variable recursion 11 

3*5*2 "Derived” Legendre functions recursion 13 

3*5*3 Recursion for terms in the radial distance 17 

3*5*4 Collection of formulas 18 

3*6 The Second Derivatives 22 

3*7 The Variational Equations 2 6 

4.0 RESULTS 29 

5*0 CONCLUSIONS 30 

APPENDIX A - GRADIENTS OF r, sin a, AND A 31 

APPENDIX B - THE CAUCHY-RIEMANN CONDITIONS 35 

APPENDIX C - ERRORS IN THE PAPER 39 

APPENDIX D - A COMPUTER PROGRAM IMPLEMENTATION 1 5 

REFERENCES U 9 


iii 



MDC W001 3 
9 February 1976 


TABLES 


Table Page 

I RECURSION SCHEMES FOR THE A 15 


n,m 


iv 



HOC WOO! 3 
9 February 1976 


FIGURES 

Figure Page 

1 Coordinates of point B . . . h 



hoc wool 3 

9 February 1976 


PINES' NONSIWULAR GRAVITATIONAL POTENTIAL 
DERIVATION, DESCRIPTION AND IMPLEMENTATION 

By Jose* L. Spencer 

McDonnell Douglas Technical Services Co., Inc. 


1.0 SUMMARY 


The possibility that the shuttle orbit er nay be required to go into polar orbits 
implies the need to derive a representation of the gravitational potential that avoids 
the usual singularity at the pole. Such a representation, including spr.eriral har- 
monics coefficients up to any order and degree, has been proposed by S. Pines in 
"Uniform Representation of the Gravitational Potential and j.cs Derivatives." 

The present note contains an engineering interpretation of and sone minor cor- 
rections to the aforementioned report by Pines. The physical meaning of the variables 
used by Pines is explained, the derivation of results is separated into smaller 
parts for easier reading, some additional recurrence relations for the "derived" 
Legendre polynomials are included and compared, and a computer program implementing 
this formulation is presented. 

\ 

Numerical experiments conducted show that the use of this representation, besides 
satisfying the requirement (removing the singularity), substantially increases the 
speed of the computation. 


2.0 INTRODUCTION 


The space shuttle is being designed with a view to its performing a multiplicity 
of tasks. Some of these may require that it be placed in a polar orbit. 

The existence of singularities in polar orbits due to the usual formulation 
of the gravitational potential makes it necessary that both ground and onboard 
software developed to support shuttle orbits use a representation of the potential 
that is free from such singularities. Pines, in reference 1, presents an alternative 
formulation that satisfies this requirement. 

The recursive algorithms proposed by Pines are stable at any order, easy to 
program, and numerically efficient. They are therefore excellent for both ground 
software (where a high order may be used) and onboard software (where an adequately 
truncated version may be obtained simply by stipulating a lower order). 

The order that must be used to meet the requirements of the onboard program 
in terms of size, speed, and accuracy is subject to study, i.r. is possible that dif- 
ferent orders may be needed for short-term and long-term s"a f * propagation. The 
ease with which the order may be changed and the variation j.1 equations obtained 
with Pines* representation of the potential make it ideal fer such a study. 


PWR0DUCIB1I ilTY OF THIt 
PAGB IB P*» 
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The formulation by Pin'. s may, then, prove to be very important for both the 
ground and the flight navigation software. For this reason, the present note has 
been written to provide (a) an engineering interpretation of reference 1, (b) a 
preliminary computer program utilizing the algorithms developed by Pines, and (c) a 
brief report on numerical results. 


3.0 DISCUSSION 


In this section, the conventional representation of the potential is presented * 
and the nature of the singularity in polar orbits is explained. 

The change of variables proposed in reference 1 is then described, and all the 
phases of the derivation are expanded; the advantages of the new representation 
are pointed out. 

It is shown how the various terms that comprise the new representation can be 
obtained recursively. For the case of the "derived" Legendre polynomials, where 
there is more than one way to obtain recurrence relations, various methods are ex- 
hibited and compared. One of these methods is recommended: Any errors that may be 

present in those terms ^rom which a new term is derived are attenuated in the process. 

Final expressions for the potential and the gravitational force are collected 
in the last subsection for easier reference. 


3.1 Statement of the Problem 

The gravitational potential of a celestial body is commonly given in spherical 
harmonics. The expression is 




cos mX 


♦ S sin ax) ] 
n,m 


( 1 ) 


where 

Li is the gravitational constant of the celestial body; 

a is the raU:us (usually the mean radius) of the celestial body; 

r is the distance between the origin of the coordinate system and the point B, 
where the potential is being evaluated; 

o is the latitude of point B; 

X is the longitude of point B; 


2 
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J are the zoual harmonics coefficients; 
n 

C and S are the coefficients of the tesseral (including the sectorial) 
n,m n,m 

harmonics; 


are the Legendre polynomials; and 

P are the associated Legendre functions. 

n,m * 

A derivation of equation (l) may be found, for instance, in reference 2. 

The choice of origin and of coordinate axes is very important in this represen- 
tation. If the origin is chosen at the celestial body’s center of mass, the constants 

J. , C and S. are all zero and the summation can start at n * 2, since 
1 l»i l *i 

the n * 1 term contributes nothing; if the z-axis is chosen to be one of the prin- 
cipal axes of inertia of the celestial body, the constants C and S 0 . are 


The angles a and A are defined as follows: a is the angle between the 
position vector of point B and its projection on the x,y plane; A is the angle 
between the positive x-axis and this projection. They are the result of the arbitrary 
choice of the x,y plane as a reference plane. 

From this arbitrary choice, there results a singularity: If point B is on the 

z-axis, A is not determined and the pctential is not defined. 

The attractive force is given by 

P = grad « = V0 - ?r + J(^l) 7{sin a) + If 7V (2) 

Here, the problem not only subsists (A is present in *|^, - ' (Via" a T* ttn ^ 

but is even aggravated by the presence of an indeterminate factor. The gradients 
(see appendix A) are 


Vr = R 


(unit vector along the position vector of point B), 

V(sin a) * — K - R 

r r 

(k is the unit vector along the positive z-axis), and 

VX = k X H 

r(l - sin 2 a) 

* * 4 

As a 90°, 1 - sin 2 a -» 0, k * R -* 0, and the factor 7 A is indeterminate . 


E £PRr,!)D 

ORIGINAL 


UUl.lTV OF TH‘- 
1VGB IS POOR 
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To obtain a representation such that the potential and the force can be unambig- 
. uously obtained for points on the z-axis, a different reference plane may be used 
(see fig. 1). 


Z 



Figure 1. - Coordinates of point B. 


But then, the problem has merely been shifted. If the y,z plane were chosen tc be 
the reference plane, the coordinates would be r, 6 * and for points cn the 
x-axis, the coordinate ^ would be undefined. If the z,x plane were chosen, the 
coordinates of the point would be r, y, and 0; and again, the 6 coordinate 
would be indeterminate for points on the y-oxis. Thus for the potential itself; to 
see what results for the force, consider the gradients of the spherical coordinates 
in the three cases mentioned. 


k 
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Plane of reference: x.y ♦ 


f * If ^ + TCsin'ai 7(sin + f? ™ 


7r = R 


7(sin a) = - k - — — a - R 
r r 


7A 


r(l - sin 2 a) 


k < R 


Plane of reference: y, 2 -*■ ? = ^ vr + 

■" — ■ ■ — • - — a** 




- v* - — — r 7(sin s) + 70/ 

3r 3(sm 6) 


( 3 ) 


7r c R 


V(sin B) = “ i - ^R ^ 
1 


Vi|> 


i x R 


Plane of reference: z,x -► 


r(l - sin 2 0) 

? *= 1^- 7r + -T7— “ — r 7(sin y) + 41* 73 
ar a(sin v) 1 36 

7r = R 


V(sin y) = - J - 2iS_ v R 
r r 




79 


r(l - sin 2 y) 


J * R 


(h) 


( 5 ) 


The problem subsists: In 7^ and 70, the sane type of iiileteminat' form that 

had been found in V\ appears again. 


3.2 The Proposed Method of Solution 

In reference 1, the method chosen to renedy this situation consists of replacing 
the spherical coordinates r, a, X (or r, £ , O' or r, > , 0 ) with another 

"distance and direction" representation of a point's position, ror this, use ic 
made of the direction cosines, which are always clearly defines for any direction 
in space. The position of each point is then given by four quantities: the three 

direction cosines that define the orientation of the position vector and r, the 
magnitude of that vector. 


V'*V (jK niI J 4 

PEVKoDU-Lm! ... ^ 

... w'-.r-' 
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Call the direction cosines s, t, and u: these are the cosines of the angles 

between the position vector and the x-axis, the y-axis, and the z-axis, respectively. 
Therefore, 

s = cos (90° - 0) * sin 0 \ 

t = cos(90° - y) = sin Y J (6) 

u = cos(90° - a) 3 sin a I 

^ A A A A A A 

On the other hand, r 3 xi + yj + zk = r(si ♦ tj + uk); so alternative ex->iesr 'on3 


C7) 


and these ere the expressions that must be used for numerical work. The attractive 
force will then be given by 

F = U Vr * otVt 7(sin 8) + aT sih r 7(sin y) * jr » u rar v(sia a) 

By using the expressions for Vr, 7(sin a), V(sin 0), and V(sin y) found in 
equations (3), (U), and (5), there results 


for s, t, and u are 


3 = * 

r 


r 


p ■ U R * juih 


sin 0 


wr(r 




or by writing s, t, and u for sin 6, sin y, and sin a, respectively, and 
collecting terras, there results 


f .fet . . tit * + + ±2±z 

\9r r 9s r 9t r 9u/ r 9s r 3t J r 9u 
which is expression (13) of reference 1. 


( 8 ) 


'The term in VX has thus been removed fro*. ?. 
yet, however; it still remains to express $ and F 


The task has not been completed 
in the new variables. 


6 
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3*3 The "Derived" Legendre Functions 

-*■ 

Certain terms in F (but not in <J>) havi Angularities that stem from another 

source: The expression for P (sin a) is, for general n and m (see, for 

n,m 

instance, refs. 3 and U), 


/•> i ,n+m 

P (sin a) = (l - sin 2 a) m — 

n ’ B 2 n n! d(sin a) n+m 


(sin 2 


if 


£ 


(1 - sin 2 a) m/2 


d®? ' in 
n 

d(sin a) 


a) 

Q 


d? n 1 (sin a) 

Bo problem exists for m > 1; but for m * 1, the first derivatives — t?— : ? — 

d(sxn a) 

are infinite at a = 190°. Write u for sin a because of equations (6) and 
differentiate. 


dP 


Oil 


du 


-u(l - u 2 ) 




,n+l 


*n , . n+1 
2 nJ du 


/, 2 K 1 /2 ,n+2 

(u 2 - l) n + ii — H -i — 1 -— ( U 2 . D 


2 n n! 


du 


n+2 


But 


du 


0+1 


(u 2 - l) n is a polynomial of degree n - 1 in u, which does not become 


n -1/2 

zero when u ® 11. The 'singularity , then, comes from the factor (l - uO ' . 


Write 


A **tn a) = 

n,m „n . ,, . xn+ro 

2 n 1 dvsin a) 


,n+m 

(sin 2 a - 1 ) D 




d(sin a) 


P (sin a) 


m n 


(9) 


and note that (1 - sin 2 a) a ^ 2 = co^ m a. Then the eumnation in in in the express. < 
of the potential becomes 


T* A (sin ci)(C 

JTi n.tn 

m-i 


cos a cos r.\ * S cos a r:n r. * } 

ii,m 


No singularities arise from differentiating A for anv value of m. (7m. A 

n,m n,m 

are polynomial".) This was done at the cost of burden. •:*; t: e terms ir. * vith the 

cos m a coefficient. The A are called by Finer* thj "cerivei" Lr g enure funct:r.n:*, 

n,m 
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3.1* The Complex Representation 

in m 

It remains, then, to express the terms cos a cos ml and cos a sin mX in 
terns of all or some of the variables r, s, t, and u. Reference 1 contains a 
very ingenious realization of this demand. 


Hote that these terms appear in a consistent form: The exponent of the cos a 

is ecpial to the coefficient of X, in all cases. A. similar behavior is to be found, 
by the application of de Moivre’s formula, in the powers of a complex number in 
polar form (see refs. 3 or U). Consider the complex number C * £ + in, where £, 

n, and c are numbers and i = / -1, and get, successively. 


in polar form, with p = ♦ n 2 and $ * 


arc tan ^ 


C ~ p(cos $ ♦ i sin ♦) 


by Euler’s formula; and 

C® ss p m (cos ♦ i sin m $) , 

by de Moivre’s formula. 

This is a new complex number w = with real and imaginary parts 

Re[ c m ] * p m cos m<& 

Imag[ ^ m ] = p m sin 


A complete analogy is found between the complex number e® and the terms 

cos®® cos mX and cos m a sin mX. These behave as the real and imaginary parts of 
the D-th power of the complex number 


C 


cos ae 


IX 


= cos a(cos X ♦ I sin X) 


( 10 ) 


Therefore, express cos a cos X and cos a sin X in terms of the variables r, 
s, t, and u (or some of then), form the complex number cos a cos X ♦ i cos "in X, 
and raise the complex number to the appropriate power to obtain the terms nee Ok* a. 


8 
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To find cos a cos A and cos o sin A in the required variables, refer back 
to figure 1 and notice the spherical triangle with sides A, a, and 90° - 6. The 
angle opposite to the side 90° - B is 90°; application of the cosine law (see, 
for instance, ref. U) gives 


cos(90° - B) * cos A cos a + sin A sin a cos 90° 


or 


sin B = cos a cos A 


( 11 ) 


Take, now, the spherical triangle with sides a, 90° - A, and 90° - y. The 
angle opposite to the side 90° - Y is 90°; apply the cosine law and find 


cos (90° - y) = cos o cos ( 90° - A) ♦ sin a sin(90° - A)cos 90° 


or * 

sin y c cos a sin A (12) 

\ 

But sin 3 and sin A are the direction cosines s and t, respectively. 

Then, s = cos a cos A, t * cos a sin A, and the complex number wanted is 

C = s + it. In reference 1, the symbols r (s,t) and i (s,t) are used to repre- 

mm 

m 

sent the real and imaginary parts, respectively, of w * £ . 


The representation of $ becomes, at this point. 


* *~jl ” £ (“) [J_ A (u) - £ A (u)(C r (s,t) + S i (s,t))]j 
r J , V r / n n,o , n,m n,m n n,m m i 

' n= l' n=1 (13) 


d“P(«) 

since, from equation (9), with u written for sin jt, A (u) = 

n ,m , n 

du 

becomes A (u) * P (u) when m * 0. 
n,o n 


There are no singularities arising from the A (u), which are the only terms 

n,m 

in u. There are none from the *ems in r. (The only one would be for r = 0, but 

the expression of the potential is only valid outside a sphere of radius a and 

center at the origin,) The terras r (s,t) and i (s,t) and their derivatives 

m r. 

with respect to s and t can be obtained from w = C . Since this is an analytic 


9 
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function of c for all the values of a considered (integer and nonnegative), the 

derivatives of all orders exist and are continuous. Differentiate v = with 
respect to C. 


dv 


a; 


tt -1 


= mlr^ts.t) ♦ ii^Cs.t)] 


But 


dv 

dc 


3s 

Si 

IS 

Bt 


+ l 


Si 
n 

9s 

3r 


a 


3t 


from the Cauchy-Riemann conditions (see appendix B), which must be satisfied because 

of the fact that w « is analytic- Since equality of complex numbers implies 
separate equalities for real and imaginary parts, the following expressions are 
obtained. 



(1U) 


These expressions give the very valuable information that i derivatives of the m-th 
terms may be obtained with no need to perform any differentiation, but only to 
multiply by m the m-lst terms. 


3.5 The Recurrence Process 

The terms in r, in u, and in (s,t) nay all be obtained recursively. The 
same is true of their derivatives of all orders. 

Recursive processes are ideally suited for use in computers. 

Tty? next subsections are devoted to the derivation of the recursion formulas 
inr the various functions involved. 

This is one of the most important parts of this note; it is upon the recursive 
i>!**’ces.. that the efficiency of the formulation depends. 


10 
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3.5*1 Complex variable recursion.- The terms required, r and i , are the 

m jjj 

real and imaginary parts of w = To obtain them recursively, form c m from 

c"- 1 . 


c = r ♦ ii 


Therefore, 


= c m-1 c 


x (r m-l + ii m-l )(s + it} 

* (s Vl - tl *.l ) + i(si m-l + tr *.l ) 


sr - - ti _ 
m-1 m-1 


i = si + tr , 
m m-1 m-1 


These are equations, (25) of reference 1. Besides these, there are two others, also 
under the number (25)* but they are wrong (see appendix C). 

The derivatives of r and i are given by equations (lU) of this note, 
m m 

To start the recurrence process, two values are needed. Since t - s + it. 


r i = s 
'1 st 

Equations (15) are valid for all integer valuer of n. The recursion could 

have been started with any other known value, such as ;° = 1, which would givt- 

r * 1, i * 0. 
o o 

The representation of the potential, equation (13), requires C R n r m ♦ 5^ 

The process of generation of the powers can be used for .this and its derivatives 

Define 


D - C r + S i 
n,n n,m m n m 


11 
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Differentiate with respect to s and to t. 


3D 3r 

— - C 


hi 


3s 


+ S 


n,m 3s a,m 3s 


m(C r . + S i . 
n,m m-1 n,m m-1 


3i 


3D 3r 

_JU£ s c + 3 5 

3t n,m 3t n,m 3t 


m(S r , - C 1 , 

q tin m— l n,m m— i 


) 


) 


Now define 


E 


n^n 


C r 4 S i 
n,m m-1 n,m m-1 


F 


n»m 


S t — C i 
n ,m m-1 n*m m— 1 


Therefore, 



mE 


n,m 


mF 


n,m 


(17) 


^Hiis Is a" 1 * that is needed to find F from it. But second derivatives are needed 
for the “ariational equations. Therefore, continue the process by differentiating 
again D with respect to s and t and defining new functions G and 

n,u» - — 

H_ then find the relations 


n.m 


n,m 



* (m - 1)G 


n,m 


= (m - 1)H 


n,m 


(18) 


AL.u, there results 


3 2 D 


3 2 D 


3s 2 


-a»s = 0 

at 2 


12 
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(D are harmonic functions) and 
n*in 


a 2 d 




n,m n,m 


a&at atas 


as vas *o be expected because of the continuity of derivatives (of all orders) of 

t. 

The functions G and H are 
n»m n*m 

G C r _ ♦ S i 0 

n*m n,m m-2 n*n m-2 


H * S r . - C i * 
n,m n,m m-2 n,3i m-2 


3*5*2 "Derived” Legendre functions recursion The A (u) are derivatives 


of the P (u). 
n 


P(U) = (u 2 - l) n 

“ 2 n! du 

by Bodrigues' formula (see refs. 3 or U) and 


n,m 


’A_ (u) = 


1 d 


n+ta 


»•“ ^n! du n+m 


(u 2 - l) n 


by definition; so 


d* 


n,m du m n 

with A n * P . 
n,0 n 

The Legendre polynomials satisfy various recurrence relations (see* for in- 
stance, ref. U). 


a) (n 4 l)P (u) - (2n ♦ l)uP (u) ♦ nP _(u) = 0 

n+1 n n-1 

b) h 1 W“» - • lr 'V ”' 1 ■ <* * i>y«> 


c) uf-lP >)] - [V -<u)] - nP <u) 
du n iu n— 1 n 


(19) 


13 
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d) 

Sr - Sr * <*» * l > p .<«> l 


(19) 

e) 

(u 2 - 1) [P (u)] = nuP (u) - nP n _(u) > 

au n n n-i 




In terms of the A (u), these relations would be wri 

tten as 



n,m 







a I ) 

(n + 1)A (u) - (2n + l)uA n (u) + nA (u) = 0 

n+1, 0 n,o 



V 

A (u) - uA 1 (u) = (n + l)A n (u) 

n+l t l n,i n,0 



c l> 

uA (u) - A .(u) * nA (u) 

Qyi n-i|i n,u 

► 

(20) 

V 

- W'*’ * (2 ° * l K,0 M 



*1> 

(uI - 1 )A d, 1 ( " ) * nuA n,0^ U * - “ 4 »-l,0 <u) 




If these relations are differentiated with respect to u, an m number of 
times, the following are obtained. 


) (n ♦ l)A .. * m(2n ♦ l)A n + (2n + l)uA + nA 


n+l,m 


n,m-l 


n,m n-1 ,m 


b-) A . * (n + m ♦ l)A + uA 

2 n+l,m+l n,m n,m+l 


c^) (n - m)A = uA - A 


n,m n,m+l n-l,m+l 


d_ ) A . , . . = (2n + 1)A + A . _ 

2 n+l,m+l n,m n-l,m+l 


e_ ) m(n - m + 1)A . = (u 2 - 1)A ♦ (2m - n)uA + nA 


n f m-l 


n,m+l 


n,m n-l,m J 


( 21 ) 


The following information is available: For m > n, A (u) = 0; for n = n, 

n ,n 

A (u) = - — — ; and for m = n - 1, A .(u) = uA (u). 
n,n 2 n nJ n,n-l n,n 

These are all straightforward consequences of the definition of A„ „ ; each A 

n ,.a n 

is a polynomial of degree n - m with only odd or only even powers of a, according 
to whether n - m is odd or even. 


Ik 
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Form a table of the A , with consideration niven all n and m from zero 
n,m — — 

on up. Let ri indicate the row and n the column of the matrix thus formed. All 
the elements of the diagonal are found by means of the expression 

A (u) = 1 * 3 * 5*7 (2n - 1 ) * " ; all elements of the upper triangular 

n,n 2 n n! 

matrix are 2ero; and the first column has the Legendre polynomials. Each element 
immediately to the left of an element of the main diagonal is the same, multiplied 
by u. 

It remains to fill the empty spaces. This can be done by any of the schemes 
that constitute the recursion relations a^ through e 2 - 

Construct the matrix, fill in the known spaces, and try to see how the various 
■\ *n - j; work. A rectangle indicates the element derived; circles indicate those it 
■ derived from (see table I). 


TABLE I. - RECURSION SCIiEftCS F(TC THE A 


! 

0 

1 

• 

2 

3 

n 

5 

6 

7 

8 

0 

( D 

0 

0 

0 

0 

0 

On 

0 

0 

1 

• 

1 

0 

0 

0 

0 

0 

0 

0 

2 

j(3u 2 -1) 

i 

a 

3 

■ 

0 

0 

0 

0 

0 . 

3 


1 

■ 

■ 

mim 

rai 

0 


0 

0 

0 

fl 

i{35u 4 -30u 2 +3) 

1 

a 


rm/mt 

8891 

105 

0 

0 

0 

0 

5 

J(63u 5 -70u 3 +15u) 





\ 945 ) 

0 

0 

0 

6 

^-(231u 6 -315u 4 +105u 2 *5) 

m 

■ 



mm 

amm 

WOT 

110395/ 

A..--- 

0 

0 

■ 

yl(429u 7 -693u 5 +315u 3 -35u) 

u 

■ 

m 


E 


135135 

0 

8 

t J h <6435u 8 -12012u 6 

+6930u 4 -1260u 2 +35) 







2027025u 

2027025 


15 
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a^) For m = 0, obtain each element of the first column in terms of the two above 
- Indicated by solid lines around the elements; for a t 0, find each element 
(not of the first column) in terms of three others: two immediately above 

and one on the previous column - indicated by dotted lines (. . . •). 

bg) Obtain elements (not of the first column) in terms of two others of the row 
above - indicated by dashed lines (- - - -). 

c^) Get elements (not of the diagonal ) in terms of two others of the following 

column - indicated by two dots-dash 

dg) Find elements (not of the first column) in terms of two others (both of previous 
rows)* one of the same column and another of the previous - indicated by dash- 
one dot . 

e^) Get elements (not of the n = 2 row) in terms of three others, Iwo of the same 

row and following columns and the other of the previous row and following 

column - indicated by one dot-two dashes ( — . — . — . — ). 

In reference 1 (where only relations and are used), relation b^ is used 
to simplify the derivatives of the potential, whereas is used for the actual re- 
cursive determination of the A • The reasoning behind this last choice is that 

n,m 

Cg is a more stable formula than b^* 


Consider the same term calculated from both formulas: "rom b„ 

A -«(n + in+l)A +uA from c_, 

n+l,m+l n,m n,m+l 2’ 


(uA 


- A 


n+l,m+l n - n n+l,ra+2 n,m+2 


,). In any error in the calculation of A 


n,n 


will be multiplied by (n + m ♦ 1) and will therefore carry a larger error into 
A . These terms are obtained from tw^ terms, one directly above and the other 

Just to the left of it; so as new terms are calculated, the values of n and m 
are increasing and the errors become progressively larger. If the term is calculated 
"by any error is divided by (n - m). Terms in the diagonil, where n - m = 0, 

are not calculated this way but serve as starting values. The recursion progresses 
to the left; and the farther away the tern is from the diagonal, the greater the 
difference between n and m and the smaller the effect of any initial error. 

For this reason , relation ^ * s to be preferred to for recursion. 

dA 

Now that the A have been found, the derivatives - - 3 must be obtained. 

n,m du 

From the definition of A , it follows that 

n,m 


r 


dA 

— = a 
du n,n+l 


( 22 ) 


1 The same argument holds without change in comparisons between and or 


a^. In e^, the effects of errors are decreased as in c^i but 
derived from three others, there are more sources of error. 


iince each term is 
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The derivative of A^ m may then be obtained by using the recursion relations, with 
no need for an actual differentiation to be performed. 


3-5-3 Recursion for terra in the radial distance .- Another quantity that may 

be obtained recursively is the term containing r. Actually, this will lead to a 
more compact representation of the potential. 

The potential was given by 


*-?!'- tffl i j , 


A (u) - 

n n,o 


m=l 


A (u)D (s,t)] 
c,m n,m 




[J A (u) 

n n,o 


m=l 


A (u)D { s ,t> ) ] 
q n ,m 


obtained by using equation (l6) in equation (13). 


The lowest D defined is D 
n,m 


1 , 1 - 


Let m = 0 in D and get 
n.m 


D _ = C r (s,t) + S i (s,t) « C because, as had been een, r = 1 ar.d 
n,o n,o o n,o o n,o o 

So if q are defined as -J^ for n = 0, the following simpler expression 
for the potential is found. 


E ♦ f; * ( tt )D is.t) 

r r W £$ n ’* n ’ m 


(23) 


A better form may be obtained by allowing n = 0 as a possible value. This 

entails defining a D =1 and A =1. In fact, A ^ =1 had already a-- 

o,o o,o o,o J 

pdared; so only the D Q Q remains to be defined. 

Since D = C , this is, in the final analysis, the only definition needed: 

C = 1. °’° 

o,o 

The expression for the potential becomes 


= Z M- 
£b rVr 


n 

Z 

m=0 


A (u)D (s,t) 
n,n n,m 


(2h) 


It is now a simple matter to find a recursive repre 


r. 

Let p * ^1 

(i)\ 

When n = 0, 

p ” when 

and 

n r 1 
it follows 

W 

that 

P„ c PP„ , • 
n n-l 

o r 


sentation for the tern in 


7 0 n-l * 


Let ~ = c , 
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As for the derivatives. 


Po 

P 1 

p 2 


*=101 


d P, 


r 

J£.« i 


dr 

dp 


o u 1 

■*5 ^ 


« »2 * 


dr 


2 u a 

r3 


ua* r 


dp 2 3ua 2 3 

* “ 4 = " a °3 

dr r 4 


PP. 


n-1 


a p n+l 


dp n _ n ♦ 1 
dr a ° n +l 


C25) 


Hote the relation 


p n = p n+l 
r a 


( 26 ) 


3.5.U Collection of formulas .- The potential can then be written as 



which is formula ( 29 ) of reference 1. 

The expression for the force was equation (7). 

+ + + 2£\r 

r 3s r 3t v r 3u \3r r 3s r }z r 3u/ 

or, in the notation of reference 1, 


V * a 2 J 


+ a-K + a, H 

0 4 


(28) 


The expressions for a^, a 2 * a^* and a^ 

tiation recursions found in subsections 3.5.1* 


c-*n now be found from the differen- 
3 * 5 * 2 , and 3*5.3. 


1 3£ 
r 3s 

Ii_ 

r 3s 


( Z p 

n=0 


n So 


*. p n 3D 

» y — y a 

£o r £o n ’ m * s 



p n+l " 

a h 


mA E 
n *m n,m 


(29) 
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a = iil 
a 2 r at 


Now 


and 


v; f mA F 

a?0 8 m=0 n ’ ran ’ n 


( 30 ) 


, a Iii 

3 r au 


Z -~£ 


JhB. v 


- a ^ * Su n*m 

n-0 m=0 


y fait f' A D 

£ 0 8 B^) n ’ m+1 


o*0 


( 31 ) 


But 


% “U- 8a l " ta 2 * ua 3 


A D = -£ ^ f* (n + 1 )A D 
dr n.m n*m a n,m n,m 

n=0 m=0 ’ n=0 m^O 


Sr 


•ua. 


= -£ fBlt y uA D 
3 n=0 8 m=0 n ’ n+1 n ’ B 


-sa. 


p n+l n 


l, - ta = - V V mA (sE ♦ tF ) 

1 2 n=0 8 g&6 n>: 3 n '° n * n 


— °n+l 


- y. — - raA [C (sr . - ti .) 
n=0 8 £=0 n * n n » B m ~ 1 . “- 1 


+ S n,m (si m-1 + tr m-l )] 

-v fail £ ^ n 

/-# a n,n n.ra 

n=0 m=0 
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So, by collecting all these, 


<Vi n 


a* *= -5] ((n + m + l)A + uA - ]D 

11 « a ft n * m n t m+l n,m 

n=0 m=0 * 


And with the help of relation b 0 from equations (2l), 


^ V* a n 

a u a " n+l,m+l n,m 

n=0 m=0 * 


(32) 


can thus be written as 


F = Y -Sii- T [mA (E i + F j) + A ,D k 
a n,m n,m n,nr n,m+l n,m 

n=o m=u 


-A ^ , ,D R] 

n+l,m+l n,m 


f 33) 


* * - 

Another way to write F, taking into account the fact that R*si+tJ+uk, is 

' 

F '(V V )i + (a 2 + a^t)J + (a^ + a J| u)k . (3 1 *) 

The main formulas are now presented together. 

Formulas for the variables: 

r = (x 2 ♦ y 2 ♦ z 2 ) 1 / 2 

•\ 

x 

s * — 
r 


*-f 


z 

u = — 
r 


(7) 
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For the recursions: 



p = OP 

n n-1 


n,n 


(2n)l 

2 n ni 


(25) 


A . = uA 
n,n-l n,n 


x ■ ■ ■ % . - _ A ) 

n,m n - m n,m+l n-l,m+l 


r o = 1 


i » 0 
o 


( 21 ) 


For the potential: 


r * 

sr 

- ti _ 

m 

ra-i 

m-1 

i = 

si n 

+ tr . 

& 

m-1 

m-1 


d = C r + s , 

n,m n,m ra-1 n,m n-1 

F * S r , - C i 

n f m n^m. n-1 n,m n-1 


D 

n t m 


C r 
n,/i n 


S i 
ri t n r 


E 

n=0 


n 

P " 

m-0 


A D 
n t n n,m 


(15) 


( 16 ) 


(27) 
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For the attractive force: 


<® o. , n 

Z-^Z 

n=0 m=0 

*» n 

Z -=r £ 

n=0 m=0 


V* p n+l 

*> n5, T 


p n+l 


mA E 

(29) 

n,m n,m 

mA F 

(30) 

n,m n,m 

A _D 

(31) 

n,m+l n,m 

^ *n+l,m+l^n,m 

• 

(32) 

* A 
(a 2 + a^t)j + (a 3 + a^ujk 

(3«») 


3.6 The Second Derivatives 

la this section, the second-order partial derivatives of $ are found. This 
can be done in several ways. In reference li the following way is chosen. 

Consider the vector 


F = Fi + FJ+Fk 
x y u z 

The "gradient of the vector" is defined to be the sum of outer products. 

^ ™ = 7(F )i^ ♦ 7(F ■ JJ 1 + 7(F )k T 

3R 2 3R * 7 2 

This is a 3 x 3 matrix, P. The vector F, in this case, was given by equation (3*0 
t - ^ ♦ a^sH + (a 2 + a^t)3 ♦ (a 3 + a^Jk 
The matrix is found to be 


3R 2 


k (a l + a u s) + V } f^ a 3 + V 5 


x~ (a, + a, s) a. ♦ a.t) |-(a- ♦ a,u) 


3y 1 


3y' 2 3y 3 


k (a i 4 v } h (a ? + V } li (a 3 + a ** u) 


pro 

11 12 13 


P P P 
21 22 r 23 


P P P 
31 *32 *33 


(35) 


and it may be seen that it is syrame+.i.c* 

22 
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T ' various elements of this matrix may be found by applying to the functions 
in parentheses the same procedure that was applied to the function The expres- 
sion for was 


V* 


dx dy dz 


the components being given by 


3x *1 b r 3s l 3r r 3s r at r 3u 


i£ = a + a t = I ii + t (M _s3i.ili.Hli) 

3y a 2 r 3t ' 3r r 3s r 3t r 3u ' 


_ a ♦ a u = — ♦ u (ii - — 12. _ £ 1£ _ H. 1±) 

3z 3 i r 3u '3r r 3s r 3t r 3u 


3 ■ da^ j 

Since — (a, + a,s) = — ^=* + sr - — + a, r^» and similarly ‘for other terms, if 

dx X ** oX oX ** dX 


will be necessary to find 


3a i 


3a, 


da ! 3% 


da, 


3a^ 3a^ 


h ’ dy 3 d2 * dx • dy * dz * dy * dz 


•, and 


1 

dz * 


(Because the matrix is symmetric, no others are needed.) As for — , J£, 

3t -du . x^y , z 

and — , recall that s » — * t = £ , and u = — and that s, t, u, and r 

are independent of one another* Therefore, 


and 


ds_ _ dt _ du ^ 1 
dx ~ dy ~ dz ~ r 


Then 


ds 9s = dt^ = 

dy dz dz 


3x * r ds S \3r 



r ds r dt r du / 


and In a similar way for the other derivatives. 
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Nov use the notation, for any function L, 


1 3L T 
r 3s * 1 


1 3L . 
r 3t * L 2 


and there results that 


1 3L r 


3L _ s 3L _ t 3Lu 3L = 

3r ’ r 3s * r r 3u ~ L U 


3a i 

-3J7 E a u + sa iu 

9a i 

W = S 12 + ta iU 


az 


s a i3 * ua H 


3*2 


W “ a 22 + ta 2i. 


aa 2 _ 

3z * a 23 + “2** 


3z * a 33 + Ua 3*» 


- — * a, , + sa, , 

3X 41 44 


3&l^ 

— = a U2 + ta w 


9a U 

- — = a. 0 + ua, | 

3Z U3 44 


Note, now, that a = a., for all i and J from 1 through k and that 
l j J * 


hi “ a 22* 


2l» 


.,..,,>11 ifx UF 

13 pooa 
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To find 


Similarly » 


5a l * a U 5 S 

p n* aT* 3 aT^ and % al and get 

p u ■ a ii * sa xu + s Ki * sa uu } * a u r 


tt u 

v + 2sa iu + s \u + r 


P 12 " P 21 - a l2 * ta iU + sa 2 U * sta U 


P 13 " P 31 * a 13 + “l>t + Sa 3U * sua UU 

p 22 * **11 + + t2a u + r 


(37) 


P„, « P-,„ * a„, ♦ ua„. + ta,. + uta 


23 32 23 2U T “ 3 b 

% 

P 33 “ a 33 + . 2Ba 3k + u2 *«* + ~ 


UU 
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The derivatives through a^ may be found by the use of formulas (17), 

( 18 ), (22), and (25). 

„ n 
p n*2 


S 11 = ~ & Z2 = , 


£ ~r £ m(m - 1)A n. m G 


n=0 ©=0 


a l2 = a oi = £ ” 2 ^ - S m ^ m - 1)A H _ m 
12 21 £b a 2 ^0 n » m n ’ ni 


n,m n,m 


a., v a - £ V nA _E 

13 31 ”- 2 Fb n » m+1 n>m 

a iu = a M = “ s r £l 

D— C 

a 23 ‘ a 32 ' £, „2 £> "*"^ 1 % 


2 ^ ^n+l Jn+l^n 

d— 0 n=0 * 


a 2U = a U2 ^ a2 


Vs n 


--£-=*£ =* 


o2 £■?„ n»= + ^ a,® 


“33 


n=0 m=0 


Su = a U3 = -£ 0 ^ S-Wm 


^n+2 . _ 

* kk * So ~ Sj a+2 ’ a+2 n - n 


i 


(38) 


3-7 The Variational Equations 

A decision on vhat order terns to include in the potential must be made after 
studying the effects on the attractive force of changes in the coefficients of the 
spherical harmonics. For th^s purpose, the partial derivatives of the force vith 
respect to those coefficients must be found* 
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For a particular coefficient C M the derivative is 

n 


3a, 


3a, 


da, 


3F ^ a i * '""2 * w ^3 * 

3C„ M * dc v M 1 + 3C. f w J * 3C„ M k * 3C, 


' R 


W N,M U,M 


"N,M 


N,M 


N,H 


( 3ft l 3a U \ " / 9a 2 * a U \ - / 9a 3 aa U \ 

» C »,M *«.») ‘ ‘ \> C N,« ' * >S,»j J * ” S,M/ 


f r i ♦ f r J ♦ f p k 

N,M,1 N,M,2 N,M,3 


(39) 


The expressions of a^, a^, a^, and a^ are given by equations (29) » (30), 

(31), and (32). 

The derivatives, then, are 


f C_ u _ a " sA N-t-l,M+i r M* 


f C , ° ~ a * tA B+l,M+i r M 


) 


f C„ . - * a (%,M+i r M “ U %+l,M+i r M 


f S a ^ MA N,M i M-l ■ sA H-*'l,M>l i M ) 

W 


f S * a M-l " tA B+l,M*l i M 

W,M,2 


y N+l 


S N,M,3 * 


^^N,M*l i M " uA N+l,M+l i M^ 


(UO) 


There are several mistakes in the expressions for these derivatives in ref- 
erence 1. For a more complete derivation of these results, see appendix C. 

Variational equations will be required for several parameters, particularly for 
the initial conditions. 

The equations of motion in an inertial system are obtained by means of a rota- 
tion matrix, from the equations relative to a rotating coordinate system. 
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Let the position vector of the point where the potential is being evaluated 

be represented by in a system of axes fixed to the celestial body and by R IW 

in an inertial system. If the rotation matrix is N, the two representations are 
connected by the equation 


NR. 


Ill 


(Ul) 


The force, which was given by = 
to the body-fixed axes), is given in the 


7 $ (where 7 means gradient with respect 

D D 

inertial system by 


-* m 

F™ * N F q 
IN B 


(U2) 


Therefore, the equations of cation of the point in the inertial system are 

<‘3> 


Consider a parameter r^. 
with respect to V ^ and to t 
that is. 


Let it be such that second-order differentiation 
can be performed in any order with the same results; 



(!*U) 


where a dot above the function stands for differentiation with respect to time. 


Then differentiate equation (1*3) with respect to 



} T 


ari f b 


and find 


+ N 


,T 

ar. 


(U5) 


T 

The matrix N depends only on the time (duration of the rotation); rep- 

resents a parameter such as the initial conditions or the coefficients of the spher- 
ical harmonics, but not the time. Therefore, 
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Then 


d / 3R I»V „T ^3 

dt V 3r i / ^ 

Carry out this last differentiation through R^. 

J V „ T i , „t. i!2nd , „T P „ is 

dt \ ar. / , s a r. ar. ar. 

\ 1/ 3R_ l l X 


(U7) 


This differentiation was carried out on the assumption that the force contains 
T t indirectly; that is, Fg contains R^, which in turn contains If F g 


should contain explicitly, a term 


f! 5 ) r , t 


must be added. 


This is- not clear in reference 1, and care must be -exercised in its use. 

The system of variational equations is, then, composed of equations (UU) and 


(■’ a, 


(U7) with, possibly, the term {N 1 ) added to equation {hi). 

explicit 

The initial conditio. 3 for this system of equations identify the values of 


— — and — at the initial time t 
dl\ dT i o 


3R. 


TH 


Solving this system of equations will give -t~ as a function of time. 

df i 


U.O RESULTS 


Pines's representation of the potential satisfies the requirement that- the 
singularities in polar orbits be removed. 

Attention is called, among the main formulas collected at th* end of section 
3.5 •*** to equation (21), which was shown to be the most stable of the recurrence 
relations for the •'derived” Legendre functions. 

A preliminary computer subroutine was written, making full us* of th*_ recursion 
relations derived by Pines. It is contained in appendix D of this note. 

James Kirkpatrick, of NASA, included this subroutine in place of the usual 
formulation of the potential in his program (which integrates the equations of 
motion with an Adams method). He noted a substantial improvement in the time of 
execution. 
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The subroutine, as shown in appendix D, makes liberal use of storage locations. 
After all the terms up to the desired order and degree have been generated and stored, 
several algebraic operations are performed upon them to obtain the components of 
the force. The terms remain in the same storage locations even after they have been 
used. 


The subroutine has since been improved with a view to saving storage locations. 
The algebraic operations are performed as soon as the necessary terns are generated, 
and only those terms that will be of some future use are retained. This has been 
the work mostly of J. Kirkpatrick. 

Although a complete evaluation of the gain in time and storage has not been 
made as yet, it is possible to say at this point that a model of the potential 
including terms up to n = m = 18 has been used with Pines f s formulation and produced 
more integration steps in less time than a model with terns up to n = m = 8 with 
the usual formulation. The integrator in which these models were included was, of 
course, the same. 


5.0 CONCLUSIONS 


The formulation proposed by Pj.nes is such that no singularities are present 
for any positions of the point. 

The recursion relations are easy to program, and the program can be made in a 
way that will result in savings in storage and time of computation; the relations 
are stable and numerically efficient at all orders. 

Generation of equations to study the effects on the attractive force of changes 
in the coefficients of the spherical harmonics or the effects on the position vector 
of changes in parameters such as those coefficients (and therefore the order of the 
potential model), initial position, or initial velocity is a simple procedure that 
makes use of the same recurrence relations. 
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APPENDIX A - GRADIENTS OF r, sin o, AND A 
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APPENDIX A 


GRADIENTS OF r, sin a, AMD A 


Consider the transformation betveen Cartesian and spherical coordinates. 


x * r cos a cos A 
y = r cos a sin A and 

z « r sin a 


r = (x 2 + y 2 + z 2 ) 1 ' 2 


z z 

sm a « - « 1 ■ ■ 1 1 ■ 

r (x 2 ♦ y 2 ♦ z 2 ) l/2 


A = arc tan ^ 

X 


The gradients of r, sin a, and A are, then. 


7r . |L t + |£ j + |r * 
3x 3y 3z 


r r r 


* ~ (x.i ♦ yj + zk) 


* — R « F 
r 


V(sin «) = «l I ♦ ikinjl * + > . ( »!" ■_ ?) “ 

3x 3y 3z 


xz 


£ _ Z5. 3 * (I _ il)K 


r r 3 


- k - — (xi + yj + zk) 
r ^ 

iJ-iiS 

r r r r 


. I i; .linjLR 

r r 
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Nov 


So 


ax ; -ax - ax r 

VXsr-l+r-J+TK 
ax ay u oz 

« — J {- *-)i f J J ♦ Ok 

2 2 X 

yi x 1 + 2 - 
1 + 2 - 2 

2 x 
V*- 


JL_i + *_ j 


x 2 ♦ y 2 x 2 + y 2 


J L i + — * ^ 


r 2 - z 2 r 2 - z 2 


-2 i + 5. 


r 2 - r 2 sin 2 a r 2 - r 2 sin 2 c 


i ♦ 


r 2 (l - sin 2 o) r^ll - sin 2 a) 


1 . 

r 2 (l - sin 2 a) 


(-yj ♦ ^d) 


k x R = 


i 

0 

x 

r 


J 

0 

2 L 

r 


k 

1 

z_ 

r 


"r i+ 7 J 


VA 


r(l - sin 2 o) 


Is. « R 
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APPENDIX B - TSE CAUCHY-RIEEAN3 CONDITIONS 
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APPENDIX B 

THE CAUCHY -R I EMANN CONDITIONS 


This derivation may be found, for instance, in reference 3. Let C = s + i 
be a complex va T able and let w « u(s,t) + iv(*,t) be a complex function of 
C t W * W(C). 

The function v(c) is said to be analytic at a point ; if the derivative 
dv • 

— exists and is finite at t and at all points of a neighborhood cf C. Tne 
rivative is defined as 


dw _ lim v(c * A;) - v ;) 
dC AC-K) A? 


where Ac = As + iAt. 

If the limit is to exist, it must be independent of the way in which A; 


dw e [u(s+As , t+At) + iv(s+As. t+At ) ] - [u(s.t) + iv(g,P] 

dC *t-K) As + iAt ' ■ 


{fu(s+As, t+At) - ufs,t)] + i[v(s+As,. t+At) - v(s,Vj}l g 
As 2 + At 2 

i{ f u( s+As , t+At) - u(s.t)] + i[v(r.+As, ^ v (s,» ) J jAt 

As 2 + At 2 


lim | 
= As-HD 
At+Oj 


Now let At * 0. Then 


dw _ lim f u( s+As ,t ) - u( s ,t ) + . v( s+As ,t ) - v^^tjl ■ >v 

dC As~*0 L As 1 As J 3s * 


Lot , now. As r 0. Then 


ii = r_,*u(s,t+At) - u ( 3 , t ) vf s , t+At ) - v(c,t) 

'P At-o L At + St 



37 


raB;PJ>BJG PAGE BLANK N<« fU-MHI 



APPENDIX B 


MOC UOOl 3 
9 February 1976 


Since equality of complex numbers requires equality of their real parts as veil 
as of their imaginary parts, the result is 



the Cauchy-Riemann conditions. 
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APPENDIX C 
ERRORS IN THE PAPER 


Although reference 1 is, in general, correct and accurate, a feu minor errors 
and misprints may lead the reader into difficulties. 

These errors and misprints are to be found in equations (10), the last two c: 
equations (25), equation (30a), and equations ( 39) • The second of equations (V>) 
is not very clear. These numbers refer to the equations in reference 1. 

Equations (10) of reference 1 read 


cos ml cos ma * r (s,t) 
m 

sin ml cos ma * i (s,t) 
m 


when they should be 


cos mX cos n a = r (s,t) 
m 

sin mX cos m a = i (s,t) 
m 


There are four equations (25), of which only the last two are wrong. Since they 
are not used afterwards in the paper, the error has no consequences. 

The functions r and i are homogeneous of degree m in s and t. Then 
o m 

Euler's theorem on homogeneous functions results in 


s 


s 


Sr Sr 

m , ^ n 
r — ♦ t tt — - el r 
3s St n 


Si Si 

m . A m . 

Ss St ” 1 n 


These could be written, with the help of the Cauchy-Eienann conditions, as 


Sr Si 

m m 

s ■— - t - — = mr 
Ss Ss m 


s 



t 


Sr 



mi 


o 


hi 
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These may be the equations that were to be presented. Instead, the last two 
of equations (25) of the reference read 

3r 3i 

8 3s * t 3t = r a 

3i a * ^a . - 
3 3s 4 * 3t * \n 

It is easy to verify that these are wrong; an exanple will suffice. Let 
, , 3r 2 

a = 2, for which r g * s z - t z , ± 2 = 2tc, = 2s, = 2s, * -2t, and 

3t 

* 2t, and obtain the results "2s 2 - 2st * s 2 - t 2 ar.i ' -2s t ♦ 2t 2 * 2s t, 
which are obviously false - 

Equation (30a) of reference 1 Just lacks a minus sign; this is corrected in 
equation ( 30b ) . 


Equations (39) of reference 1 have several mistakes. As shown in section 3*7 
of this note, the derivatives are 




3a 1 3a l» 


3a i 3a b 


3C ^ ® * f. * jc * ® jc 

B.M S,M H,M,1 ",M H,M 


3a 2 3a U 

3C„ .. 4 * 30 


3a^ 3a^ 

asT 7. 4 z 3s“ 


M, ,2 N,M H,M H,M,2 S,H 


3a, 


3C„ u 3C. 


3a u 


3a 


3 3a U 

♦ u 


n,m,3 n,m ° nj* n,m,3 ° :;,m n,m 


3S.. * 3S m 


These derivatives are to be obtained from the expressions of a^, a^, a^, 

and a^; namely, 

v -Sii y aA (C r ,♦ S i ,) 

1 ^ a n,m n,m n-1 a— 1 

n-0 m=0 

- P_., n 

a « T -2i± Y mA (S r , - C i J 

2 a n *m n,m n-1 n,r. a-1 

n-0 m— J 


°n+l n 


« - y - 2 ^ y A ,(C r S i ) 

* &> m £s& n ’ n+1 n,a n,a a 


e n-»l n 


a. . -V -Sli y A , Ac r ♦ S i ) 

h a n+l,ra+l n,n s n,m m 

n=0 m=0 


J.2 
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Consider a specific pair of coefficients C„ and S„ ... They appeal* only 
once in the expressions of a^, and a^. The term may be written 
so as to show separately the part containing these coefficients. 


"- 1 Vi A 


«W M - x 


a = V -2I±. V cA E + y nsA E 

^ a n ’ m n ’ B a m^o* u ’ m 




N+l 


* a * .a 


Pm* i 5 a> 0-4 i n 

_NHja y . E + y -BlL ymA E 

a *,m • 4i, a ^ n,n n*n 

m=M+l * * n=N+l m=Q 

Similarly for a 2 » a^» and (but showing only the pertinent parts). 


a 2 = — * a 


»♦! 


y N+l 




TM 


N*1 


a 3 = ••• * a A BJM C I.R r M * a ^1 * 


K+l 


N+l 


~ a ,K^1 S N * 


vlth the result that 

9a l P N+1 


3C W „ ’ a ’ 

W 


9a , 


ac 


N,M 


*M+1 lt _ 


**3 °N+1 


ac., 


9a 


as 


1 - p bHma « 


N,H 

0 *'N+1 


3S h,m 

S ■ ■ ■ ■ 

a 

9a 3 

P N+1 

ss • — 




9 a. 


ac, 


N,M 


JHa r ; 

a ti+l,M+l M 


9a, 


°W+1 


as. 


N,M 


A s+1,M+1 1 M 


Knowing these, it is an easy matter to arrive at the correct equations, which 
may be found in section 3*7 of this note as equations (^0). 

The lack of clarity of the second of equations (U 5 ) of reference 1 has been 
discussed at the end of section 3.7 of the present note. 
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APPENDIX D 

A COMPUTER PROGRAM IMPLEMENTATION 


SUBROUTINE GRAVPOT (NAXO,X,Y,Z,FX,FY,FZ) 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DIMENSION CREAL(flAXO), CIMAG(NAXO), RHO(NAXO), 

1 A(HAXO.NAXO) ,D(riAXO,NAXO) ,E(NAXO,NAXO) ,F(KAXO,NAXO) 

COMMON/ CONST / C ( N AXO ,N AXO ) ( NAXO ,H AXO ) .ZONAL(NAXO) , 

1 TERRAD, TMU 
C 

C NAXO IS MAXO+1 , WHERE MAXO IS THE ORDER TO WHICH THE 
C POTENTIAL IS DESIRED. 

C 

C TERRAD IS THE EARTH’S RADIUS, MU THE GRAVITATIONAL CONSTANT, 
C C(N,M), S(N,M) ANO ZONAL(M) THE COEFFICIENTS OF THE 
C SPHERICAL HARMONICS, A(N,M) THE DERIVED LEGENORF 
C FUNCTIONS, X, Y ANO Z THE COORDINATES OF THE POINT, 

C D(N.M), E(N,M) AND F(N,M) THE FUNCTIONS OF THE 
C COMPLEX VARIABLE ZET A= CREAL+ 1 *C I MAG AND 
C RHO IS THE FUNCTION OF THE RADIAL DISTANCE. FX, 

C FY ANO FZ ARE THE COMPONENTS OF THE FORCE. 

C 

C GET THE DIRECTION COSINES ESS, T AND U. 

C 

MAXO=NAXO-1 

RINV=1 . ODG/SQRT ( X*X+Y*Y +Z*Z ) 

ESS=X*RINV 

T=Y*RINV 

U*Z*RINV 

C 

C GENERATE THE FUNCTIONS CREAL, CIMAG, A, D, E, F AND RHO. 

C 

RO=TERRAD*RINV 
RHOZERO s TM"*RHOZERO 
RHO(1 )=RO*RHOZERO 
CREAL(1 )=ESS 
CIMAG(1 )=T 
D(1 ,1 )=O.ODO 
E(1 ,1 )=O.ODO 
F(1 ,1 )=O.ODO 
A(1 ,1 )=1 .ODO 
DO 1 1*2, NAXO 

RHO(I)=RO+RHO(I-l) 


*»7 


FAG* 


Vt/.NX 


hij*.:.. 
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CREAl ( I )=ESS*CREAL ( I - 1 ) -T*CIMAG{ I - 1 ) 

CIMAG( 1 )=E$S*CIMAG(I-1 )+T*CREAL( 1-1 ) 
D(I,1)=ESS*C(I,1)+T*S(I.1) 

E(1 ,1)=C(I,1) 

F(I.1)-S(M) 

A( I ,1)* (2*1*1 )*A( I-l.I-l ) 

A(I,I-1)=U*A(I,I ) 

00 1 K*2,I 

IF (I.EQ.NAXO) GO TO 9 

D( I ,K)=C(I ,K)*CREAL(K)*S(I ,K)*CIMAG(K) 

E(I,K)=C(I,K)*CREAL{K-1)+S(I,K)*CIMAG(K-1) 

F( I ,K)=S(I ,K)*CREAL(K-1 )-C( I ,K)*CIMAG(K-1 ) 

9 CONTINUE 

IF <I.EQ.2) GO TO 10 
L=I-2 

DO 1 J=1,L 

A(I,I-.M)=(U*A(I,I.J)-A(I-1,I-J))/(J+1) 

1 COfPINUE 
10 CONTINUE 

C NO!) COWUTE THE AUXILIARY QUANTITIES Al, A2, A3 AND 

C A4, NEEDED TO FIND THE COMPONENTS OF THE FORCE. 

A1=0.0D0 

A2-O.ODO 

A3=O.ODO 
A4=RH0ZER0*RINV 
DO 2 N=2,MAX0 
FAC 1=0. 000 
FAC2=0.000 
FAC3=A(N,1)*Z0NAL(N) 

FAC4=A(H+1 ,1 )*ZOtlAL(;N) 

00 3 M=1 ,H 

FAC1 =F AC1 +M* A( N ,M ) *£ ( N ,M ) 

FAC2= F AC2+M* A ( N ,M } *F ( N ,M) 

FAC3=FAC3*A(N,MO )*D(N,M) 

3 FAC4=FAC4+A(fl+1 ,M+1 )*D(N,M) 

A1=A1+RINV*RH0(N)*FAC1 
A2=A2+RIMV*RH0(N)*FAC2 
A3= A3+R INV* RHO ( N ) *F AC3 
A4=A4+RINV* RHO (N )*FAC4 

2 CONTINUE 
FX=A1-ESS*A4 
FY=A2-T*A4 
FZ=A3-U*A4 
RETURN 

END 
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